ROOT := "/home/kohel/teaching/2002.MATH3024/Solutions/";
LIBS := ROOT * "digraph_frequencies/";
TXTS := "/home/kohel/literature/PoeEdgarAllan/";
AZ := [ CodeToString(64+i) : i in [1..26] ];
A2 := [ X*Y : X, Y in AZ ];
AA := [ X*X : X in AZ ];
S := StripEncoding(Read(TXTS * "blackcat.txt")); 
F1 := FrequencyDistribution(S);
F2 := DigraphFrequencyDistribution(S);
function SprintReal(r,n)
    error if r gt 1, "Argument 1 must be at most 1";
    rs := Sprint(r);
    if &or[ "e" eq rs[i] : i in [1..#rs] ] then
	m := Ceiling(-Log(10,r))-1;
	ss := Sprint(r*10^m);
	rs := "0." * &cat[ "0" : i in [1..m] ];
	rs *:= &*[ Strings() | ss[i] : i in [3..#ss] ];
    end if;
    if rs eq "0" then
	rs := "0." * &*[ "0" : i in [1..n] ];
    elif #rs gt n+2 then
	rs := &*[ rs[i] : i in [1..n+2] ];
    end if;
    while #rs lt n+2 do rs *:= "0"; end while;
    return rs;
end function;
function SprintVector(v,n)
    vs := "(" * SprintReal(v[1],n);
    for i in [2..Degree(Parent(v))] do
	vs *:= ",\ " * SprintReal(v[i],n);
    end for;
    return vs * ")";
end function;
MS := Strings(); 
V := [ Vector([ RR | F1[i],
       F1[i] eq 0 select 0 else F2[i,i]/F1[i] ]) : i in [1..26] ];
eps0 := RealField(8)!0;
eps1 := RealField(8)!0.00001;
eps2 := RealField(8)!0.00003;
eps3 := RealField(8)!0.00010;
eps4 := RealField(8)!0.00032;
N := [ < AA[i], V[i], 
     {@ AA[j] : j in [1..26] | 
         Norm(V[i]-V[j]) ge eps0 and 
         Norm(V[i]-V[j]) lt eps1 @},
     {@ AA[j] : j in [1..#V] | 
         Norm(V[i]-V[j]) ge eps1 and 
         Norm(V[i]-V[j]) lt eps2 @},
     {@ AA[j] : j in [1..#V] | 
	 Norm(V[i]-V[j]) ge eps2 and 
	 Norm(V[i]-V[j]) lt eps3 @},
     {@ AA[j] : j in [1..#V] | 
	 Norm(V[i]-V[j]) ge eps3 and 
	 Norm(V[i]-V[j]) lt eps4 @}>
   : i in [1..26] ];
/* */
FILE := Open(LIBS * "digraph_homeo_freqs.tex","w"); 
Puts(FILE,"$$");
Puts(FILE,"\\begin{array}{ccccccc}");
Puts(FILE,"X & (P(X),P(XX|X)) & " *
    SprintReal(eps1,5) * " &" * SprintReal(eps2,5) * " &" * 
    SprintReal(eps3,5) * " &" * SprintReal(eps4,5) * "\\\\ \\hline");
for X in N do
    line := "{\\tt " * X[1][1] * "} & " * SprintVector(X[2],6) * " & ";
    for k in [3..#X-1] do
	if #X[k] ne 0 then
	    line *:= "{\\tt " *
	        &*[ MS | X[k][i][1] * "," : i in [1..#X[k]-1] ] *
	        X[k][#X[k]][1] * "} & ";
	else
	    line *:= " & "; 
	end if;
    end for;
    if #X[#X] ne 0 then
	line *:= "{\\tt " *
	    &*[ MS | X[#X][i][1] * "," : i in [1..#X[#X]-1] ] *
	    X[#X][#X[#X]][1] * "}\\\\";
    else
	line *:= " \\\\ "; 
    end if;
    Puts(FILE,line);
end for;
Puts(FILE,"\\end{array}");
Puts(FILE,"$$");
delete FILE;
/* */
// Projection to 2-dimensional graph: (x,y,z) -> (u,v), where
//     u := (y sin(theta) + z cos(theta)) sin(phi),
//     v := x cos(phi) - y cos(theta) + z sin(theta).
phi := Pi(RealField())/12;
theta := Pi(RealField())/6;
s1 := RealField(8)!Cos(phi); t1 := RealField(8)!Cos(theta); 
s2 := RealField(8)!Sin(phi); t2 := RealField(8)!Sin(theta); 
raxis := 0.19;
xaxis := Sprintf("(%o,%o)", "0.000000", SprintReal(raxis*s1,6));
yaxis := Sprintf("(%o,%o)",
    SprintReal(raxis*t2,6), SprintReal(-raxis*s2*t1,6));
zaxis := Sprintf("(%o,%o)",
    SprintReal(raxis*t1,6), SprintReal(raxis*s2*t2,6));
// for k in [1,5,8,9,12,13,14,15,18,19,20,21] do
for k in [1..26] do
    C := CodeToString(64+k);
    VC := [ Vector([RR|F1[i],F2[i,k]/F1[k],F2[k,i]/F1[k]]) : i in [1..26] ];
    eps0 := RealField(8)!0;
    eps1 := RealField(8)!0.00001;
    eps2 := RealField(8)!0.00005;
    eps3 := RealField(8)!0.00010;
    eps4 := RealField(8)!0.00015;
    NC := [ < A2[k+(i-1)*26], VC[i], 
       {@ A2[k+(j-1)*26] : j in [1..26] | 
           delta ge eps0 and delta lt eps1 where
           delta := Norm(VC[i]-VC[j]) @},
       {@ A2[k+(j-1)*26] : j in [1..26] | 
           delta ge eps1 and delta lt eps2 where
           delta := Norm(VC[i]-VC[j]) @},
       {@ A2[k+(j-1)*26] : j in [1..26] | 
           delta ge eps2 and delta lt eps3 where
           delta := Norm(VC[i]-VC[j]) @},
       {@ A2[k+(j-1)*26] : j in [1..26] | 
           delta ge eps3 and delta lt eps4 where
           delta := Norm(VC[i]-VC[j]) @}>
        : i in [1..26] ];
    FILE := Open(LIBS * "digraph_mono" * C * "_freqs.tex","w"); 
    Puts(FILE,"$$");
    Puts(FILE,"\\begin{array}{ccccccc}");
    Puts(FILE,"X & " * 
        "(P(X),P(X{\\tt " * C * "}|{\\tt " * C * "}))," * 
        "P({\\tt " * C * "}X|{\\tt " * C * "})) & " *
//        SprintReal(eps1,5) * " &" * SprintReal(eps2,5) * " &" * 
//        SprintReal(eps3,5) * " &" * SprintReal(eps4,5) * "\\\\ \\hline");
          " 1\\epsilon &  5\\epsilon & " *
          "10\\epsilon & 15\\epsilon \\\\ \\hline");
    for X in NC do
        line := "{\\tt " * X[1][1] * "} & " * SprintVector(X[2],5) * " & ";
        for k in [3..#X-1] do
	    if #X[k] ne 0 then
	        line *:= "{\\tt " *
	            &*[ MS | X[k][i][1] * "," : i in [1..#X[k]-1] ] *
  	            X[k][#X[k]][1] * "} & ";
            else
	        line *:= " & "; 
   	    end if;
        end for;
        if #X[#X] ne 0 then
	    line *:= "{\\tt " *
	        &*[ MS | X[#X][i][1] * "," : i in [1..#X[#X]-1] ] *
                X[#X][#X[#X]][1] * "}\\\\";
        else
	    line *:= " \\\\ "; 
        end if;
        Puts(FILE,line);
    end for;
    Puts(FILE,"\\end{array}");
    Puts(FILE,"$$");
    delete FILE;
    /* */
    FILE := Open(LIBS * "digraph_mono" * C * "_graph.tex","w"); 
    Puts(FILE,"$$");
    Puts(FILE,"\\begin{xy}");
    Puts(FILE,"0;<24cm,0cm>:<0cm,24cm>::");
    Puts(FILE,"@i;");
    Puts(FILE,"(0,0)=*={o} +\"o\";");
    Puts(FILE,xaxis * "=*={x} +\"x\" *{x} **@{-} ?>*@{>};");
    Puts(FILE,"(0,0)=*={o} +\"o\";");
    Puts(FILE,yaxis * "=*={y} +\"y\" *{y} **@{-} ?>*@{>};");
    Puts(FILE,"@i;");
    Puts(FILE,"(0,0)=*={o} +\"o\";");
    Puts(FILE,zaxis * "=*={z} +\"z\" *{z} **@{-} ?>*@{>};");
    for X in NC do
        C2 := X[1][1];
        C2tt := "{\\tt " * C2 * "}";
        Puts(FILE,"@i;");
        Puts(FILE,"(0,0)=*={o} +\"o\";");
        Puts(FILE,Sprintf("(%o,%o)=*={%o} +\"%o\" *{%o} **@{-} ?>*@{>};",
    	    SprintReal(t2*X[2][2]+t1*X[2][3],6),
	    SprintReal(s1*X[2][1]+s2*(-t1*X[2][2]+t2*X[2][3]),6),
            C2,C2,C2tt));
     end for;
     Puts(FILE,"\\end{xy}");
     Puts(FILE,"$$");
     delete FILE;
end for;
/* */
FILE := Open(LIBS * "digraph_homeo_graph.tex","w"); 
Puts(FILE,"$$");
Puts(FILE,"\\begin{xy}");
Puts(FILE,"0;<36cm,0cm>:<0cm,36cm>::");
Puts(FILE,"(0,0)=*={o} +\"o\";");
Puts(FILE,"(0.13,0.00)=*={x} +\"x\" *{x} **@{-} ?>*@{>};");
Puts(FILE,"@i;");
Puts(FILE,"(0,0)=*={o} +\"o\";");
Puts(FILE,"(0.00,0.14)=*={y} +\"y\" *{y} **@{-} ?>*@{>};");
for X in N do
    C := X[1][1];
    Ctt := "{\\tt " * C * "}";
    Puts(FILE,"@i;");
    Puts(FILE,"(0,0)=*={o} +\"o\";");
    Puts(FILE,Sprintf("(%o,%o)=*={%o} +\"%o\" *{%o} **@{-} ?>*@{>};",
        SprintReal(X[2][1],6),SprintReal(X[2][2],6),C,C,Ctt));
end for;
Puts(FILE,"\\end{xy}");
Puts(FILE,"$$");
delete FILE;
/* */
function ProbabilityVector(F1,F2,XY);
    i := StringToCode(XY[1])-64;
    j := StringToCode(XY[2])-64;
    return Vector([ RealField(8) |
    F1[i], F1[j], F2[i,j], 
    F1[i] eq 0 select 0 else F2[i,j]/F1[i],
    F1[j] eq 0 select 0 else F2[i,j]/F1[j],
    F1[i] eq 0 select 0 else F2[j,i]/F1[i],
    F1[j] eq 0 select 0 else F2[j,i]/F1[j]]);
end function;
FILE := Open(LIBS * "digraph_hetro_freqs.tex","w"); 
for Isol in [ Split(Read(LIBS * "the_list_sub1")),
              Split(Read(LIBS * "the_list_sub2")),
              Split(Read(LIBS * "the_list_sub3")),
              Split(Read(LIBS * "the_list_sub4")) ] do
    Puts(FILE,"$$");
    Puts(FILE,"\\begin{array}{cc}");
    Puts(FILE,"XY & (P(X),P(Y),P(XY),P(XY|X),P(XY|Y),P(YX|Y),P(YX|X))\\\\\\\hline");
    for XY in Isol do
	VS :=SprintVector(ProbabilityVector(F1,F2,XY),6); 
	Puts(FILE,"{\\tt " * XY * "}" * " & " * VS * "\\\\");
    end for;
    Puts(FILE,"\\end{array}");
    Puts(FILE,"$$");
end for;
delete FILE;
/* */
FILE := Open(LIBS * "digraph_hetro_ngbrs.tex","w"); 
Isol := Split(Read(LIBS * "the_list_sub1"));
VV := [ ProbabilityVector(F1,F2,XY) : XY in A2 ];
eps0 := RealField(8)!0;
eps1 := RealField(8)!0.001;
eps2 := RealField(8)!0.005;
eps3 := RealField(8)!0.010;
eps4 := RealField(8)!0.015;
eps5 := RealField(8)!0.020;
NN := [ < XY, 
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps0 and Norm(VV[i]-VV[j]) lt eps1 @},
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps1 and Norm(VV[i]-VV[j]) lt eps2 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps2 and Norm(VV[i]-VV[j]) lt eps3 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps3 and Norm(VV[i]-VV[j]) lt eps4 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps4 and Norm(VV[i]-VV[j]) lt eps5 @} >
     where i := Index(AZ,XY[2])+26*(Index(AZ,XY[1])-1) : XY in Isol ];
Puts(FILE,"$$");
Puts(FILE,"\\begin{array}{ccccccc}");
Puts(FILE,"XY &\n" * 
    "\\hspace{0.3cm}" * SprintReal(eps1,4) * "\\hspace{0.3cm} &\n" * 
    "\\hspace{0.0cm}" * SprintReal(eps2,4) * "\\hspace{0.0cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps3,4) * "\\hspace{0.6cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps4,4) * "\\hspace{0.6cm} &\n " * 
    "\\hspace{0.7cm}" * SprintReal(eps5,4) * "\\hspace{0.7cm}\\\\ \\hline");
for X in NN do
    line := "{\\tt " * X[1] * "}";
    for k in [2..#X] do
	ss := Sprint(X[k]);
	line *:= " & {\\tt " * &*[ ss[i] : i in [3..#ss-2] ] * "}";
    end for;
    line *:= "\\\\";
    Puts(FILE,line);
end for;
Puts(FILE,"\\end{array}");
Puts(FILE,"$$");
/* */
Isol := Split(Read(LIBS * "the_list_sub2"));
eps0 := RealField(8)!0;
eps1 := RealField(8)!0.0010;
eps2 := RealField(8)!0.0025;
eps3 := RealField(8)!0.0050;
eps4 := RealField(8)!0.0075;
eps5 := RealField(8)!0.0100;
NN := [ < XY, 
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps0 and Norm(VV[i]-VV[j]) lt eps1 @},
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps1 and Norm(VV[i]-VV[j]) lt eps2 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps2 and Norm(VV[i]-VV[j]) lt eps3 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps3 and Norm(VV[i]-VV[j]) lt eps4 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps4 and Norm(VV[i]-VV[j]) lt eps5 @} >
     where i := Index(AZ,XY[2])+26*(Index(AZ,XY[1])-1) : XY in Isol ];
Puts(FILE,"$$");
Puts(FILE,"\\begin{array}{ccccccc}");
Puts(FILE,"XY &\n" * 
    "\\hspace{0.3cm}" * SprintReal(eps1,4) * "\\hspace{0.3cm} &\n" * 
    "\\hspace{0.0cm}" * SprintReal(eps2,4) * "\\hspace{0.0cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps3,4) * "\\hspace{0.6cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps4,4) * "\\hspace{0.6cm} &\n " * 
    "\\hspace{0.7cm}" * SprintReal(eps5,4) * "\\hspace{0.7cm}\\\\ \\hline");
for X in NN do
    line := "{\\tt " * X[1] * "}";
    for k in [2..#X] do
	ss := Sprint(X[k]);
	line *:= " & {\\tt " * &*[ ss[i] : i in [3..#ss-2] ] * "}";
    end for;
    line *:= "\\\\";
    Puts(FILE,line);
end for;
Puts(FILE,"\\end{array}");
Puts(FILE,"$$");
/* */
Isol := Split(Read(LIBS * "the_list_sub3"));
eps0 := RealField(8)!0;
eps1 := RealField(8)!0.0010;
eps2 := RealField(8)!0.0020;
eps3 := RealField(8)!0.0030;
eps4 := RealField(8)!0.0040;
eps5 := RealField(8)!0.0050;
NN := [ < XY, 
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps0 and Norm(VV[i]-VV[j]) lt eps1 @},
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps1 and Norm(VV[i]-VV[j]) lt eps2 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps2 and Norm(VV[i]-VV[j]) lt eps3 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps3 and Norm(VV[i]-VV[j]) lt eps4 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps4 and Norm(VV[i]-VV[j]) lt eps5 @} >
     where i := Index(AZ,XY[2])+26*(Index(AZ,XY[1])-1) : XY in Isol ];
Puts(FILE,"$$");
Puts(FILE,"\\begin{array}{ccccccc}");
Puts(FILE,"XY &\n" * 
    "\\hspace{0.3cm}" * SprintReal(eps1,4) * "\\hspace{0.3cm} &\n" * 
    "\\hspace{0.0cm}" * SprintReal(eps2,4) * "\\hspace{0.0cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps3,4) * "\\hspace{0.6cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps4,4) * "\\hspace{0.6cm} &\n " * 
    "\\hspace{0.7cm}" * SprintReal(eps5,4) * "\\hspace{0.7cm}\\\\ \\hline");
for X in NN do
    line := "{\\tt " * X[1] * "}";
    for k in [2..#X] do
	ss := Sprint(X[k]);
	line *:= " & {\\tt " * &*[ ss[i] : i in [3..#ss-2] ] * "}";
    end for;
    line *:= "\\\\";
    Puts(FILE,line);
end for;
Puts(FILE,"\\end{array}");
Puts(FILE,"$$");
/* */
Isol := Split(Read(LIBS * "the_list_sub4"));
eps0 := RealField(8)!0;
eps1 := RealField(8)!0.0010;
eps2 := RealField(8)!0.0015;
eps3 := RealField(8)!0.0020;
eps4 := RealField(8)!0.0025;
eps5 := RealField(8)!0.0030;
NN := [ < XY, 
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps0 and Norm(VV[i]-VV[j]) lt eps1 @},
     {@ A2[j] : j in [1..#VV] | 
         Norm(VV[i]-VV[j]) ge eps1 and Norm(VV[i]-VV[j]) lt eps2 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps2 and Norm(VV[i]-VV[j]) lt eps3 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps3 and Norm(VV[i]-VV[j]) lt eps4 @},
     {@ A2[j] : j in [1..#VV] | 
	 Norm(VV[i]-VV[j]) ge eps4 and Norm(VV[i]-VV[j]) lt eps5 @} >
     where i := Index(AZ,XY[2])+26*(Index(AZ,XY[1])-1) : XY in Isol ];
Puts(FILE,"$$");
Puts(FILE,"\\begin{array}{ccccccc}");
Puts(FILE,"XY &\n" * 
    "\\hspace{0.3cm}" * SprintReal(eps1,4) * "\\hspace{0.3cm} &\n" * 
    "\\hspace{0.0cm}" * SprintReal(eps2,4) * "\\hspace{0.0cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps3,4) * "\\hspace{0.6cm} &\n" * 
    "\\hspace{0.6cm}" * SprintReal(eps4,4) * "\\hspace{0.6cm} &\n " * 
    "\\hspace{0.7cm}" * SprintReal(eps5,4) * "\\hspace{0.7cm}\\\\ \\hline");
for X in NN do
    line := "{\\tt " * X[1] * "}";
    for k in [2..#X] do
	ss := Sprint(X[k]);
	line *:= " & {\\tt " * &*[ ss[i] : i in [3..#ss-2] ] * "}";
    end for;
    line *:= "\\\\";
    Puts(FILE,line);
end for;
Puts(FILE,"\\end{array}");
Puts(FILE,"$$");
delete FILE;
